import os
import requests
import boto3
from osgeo import gdal
import rasterio as rio
from rasterio.session import AWSSession
import rioxarray
import hvplot.xarray
import holoviews as hvAccessing Cloud Optimized GeoTIFF (COG) from Earthdata Cloud
Summary
Harmonized Landsat Sentinel-2 (HLS) Operational Land Imager Surface Reflectance and TOA Brightness Daily Global 30m v2.0 (L30) (10.5067/HLS/HLSL30.002)
Requirements
AWS intance running in us-west 2
Earthdata Login
.netrc file
Learning Objectives
- S3 Direct Access
Cloud Optimized GeoTIFF (COG)
Using Harmonized Landsat Sentinel-2 (HLS) version 2.0
Import Packages
Single File S3 Direct Access
Temporary Credentials
s3_cred_endpoint = {
'podaac':'https://archive.podaac.earthdata.nasa.gov/s3credentials',
'gesdisc': 'https://data.gesdisc.earthdata.nasa.gov/s3credentials',
'lpdaac':'https://data.lpdaac.earthdatacloud.nasa.gov/s3credentials',
'ornldaac': 'https://data.ornldaac.earthdata.nasa.gov/s3credentials',
'ghrcdaac': 'https://data.ghrc.earthdata.nasa.gov/s3credentials'
}def get_temp_creds(provider):
return requests.get(s3_cred_endpoint[provider]).json()temp_creds_req = get_temp_creds('lpdaac')
#temp_creds_reqsession = boto3.Session(aws_access_key_id=temp_creds_req['accessKeyId'],
aws_secret_access_key=temp_creds_req['secretAccessKey'],
aws_session_token=temp_creds_req['sessionToken'],
region_name='us-west-2')GDAL Configurations
GDAL is a foundational piece of geospatial software that is leveraged by several popular open-source, and closed, geospatial software. The rasterio package is no exception. Rasterio leverages GDAL to, among other things, read and write raster data files, e.g., GeoTIFFs/Cloud Optimized GeoTIFFs. To read remote files, i.e., files/objects stored in the cloud, GDAL uses its Virtual File System API. In a perfect world, one would be able to point a Virtual File System (there are several) at a remote data asset and have the asset retrieved, but that is not always the case. GDAL has a host of configurations/environmental variables that adjust its behavior to, for example, make a request more performant or to pass AWS credentials to the distribution system. Below, we’ll identify the evironmental variables that will help us get our data from cloud
rio_env = rio.Env(AWSSession(session),
GDAL_DISABLE_READDIR_ON_OPEN='TRUE',
GDAL_HTTP_COOKIEFILE=os.path.expanduser('~/cookies.txt'),
GDAL_HTTP_COOKIEJAR=os.path.expanduser('~/cookies.txt'))
rio_env.__enter__()<rasterio.env.Env at 0x7f7ad2881d60>
s3_url = 's3://lp-prod-protected/HLSL30.020/HLS.L30.T11SQA.2021333T181532.v2.0/HLS.L30.T11SQA.2021333T181532.v2.0.B04.tif'da = rioxarray.open_rasterio(s3_url)da<xarray.DataArray (band: 1, y: 3660, x: 3660)>
[13395600 values with dtype=int16]
Coordinates:
* band (band) int64 1
* x (x) float64 7e+05 7e+05 7e+05 ... 8.097e+05 8.097e+05 8.097e+05
* y (y) float64 4.1e+06 4.1e+06 4.1e+06 ... 3.99e+06 3.99e+06
spatial_ref int64 0
Attributes:
_FillValue: -9999.0
scale_factor: 0.0001
add_offset: 0.0
long_name: Red- band: 1
- y: 3660
- x: 3660
- ...
[13395600 values with dtype=int16]
- band(band)int641
array([1])
- x(x)float647e+05 7e+05 ... 8.097e+05 8.097e+05
array([699975., 700005., 700035., ..., 809685., 809715., 809745.])
- y(y)float644.1e+06 4.1e+06 ... 3.99e+06
array([4100025., 4099995., 4099965., ..., 3990315., 3990285., 3990255.])
- spatial_ref()int640
- crs_wkt :
- PROJCS["UTM Zone 11, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
- semi_major_axis :
- 6378137.0
- semi_minor_axis :
- 6356752.314245179
- inverse_flattening :
- 298.257223563
- reference_ellipsoid_name :
- WGS 84
- longitude_of_prime_meridian :
- 0.0
- prime_meridian_name :
- Greenwich
- geographic_crs_name :
- Unknown datum based upon the WGS 84 ellipsoid
- horizontal_datum_name :
- Not_specified_based_on_WGS_84_spheroid
- projected_crs_name :
- UTM Zone 11, Northern Hemisphere
- grid_mapping_name :
- transverse_mercator
- latitude_of_projection_origin :
- 0.0
- longitude_of_central_meridian :
- -117.0
- false_easting :
- 500000.0
- false_northing :
- 0.0
- scale_factor_at_central_meridian :
- 0.9996
- spatial_ref :
- PROJCS["UTM Zone 11, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
- GeoTransform :
- 699960.0 30.0 0.0 4100040.0 0.0 -30.0
array(0)
- _FillValue :
- -9999.0
- scale_factor :
- 0.0001
- add_offset :
- 0.0
- long_name :
- Red
da_red = da.squeeze('band', drop=True)
da_red<xarray.DataArray (y: 3660, x: 3660)>
[13395600 values with dtype=int16]
Coordinates:
* x (x) float64 7e+05 7e+05 7e+05 ... 8.097e+05 8.097e+05 8.097e+05
* y (y) float64 4.1e+06 4.1e+06 4.1e+06 ... 3.99e+06 3.99e+06
spatial_ref int64 0
Attributes:
_FillValue: -9999.0
scale_factor: 0.0001
add_offset: 0.0
long_name: Red- y: 3660
- x: 3660
- ...
[13395600 values with dtype=int16]
- x(x)float647e+05 7e+05 ... 8.097e+05 8.097e+05
array([699975., 700005., 700035., ..., 809685., 809715., 809745.])
- y(y)float644.1e+06 4.1e+06 ... 3.99e+06
array([4100025., 4099995., 4099965., ..., 3990315., 3990285., 3990255.])
- spatial_ref()int640
- crs_wkt :
- PROJCS["UTM Zone 11, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
- semi_major_axis :
- 6378137.0
- semi_minor_axis :
- 6356752.314245179
- inverse_flattening :
- 298.257223563
- reference_ellipsoid_name :
- WGS 84
- longitude_of_prime_meridian :
- 0.0
- prime_meridian_name :
- Greenwich
- geographic_crs_name :
- Unknown datum based upon the WGS 84 ellipsoid
- horizontal_datum_name :
- Not_specified_based_on_WGS_84_spheroid
- projected_crs_name :
- UTM Zone 11, Northern Hemisphere
- grid_mapping_name :
- transverse_mercator
- latitude_of_projection_origin :
- 0.0
- longitude_of_central_meridian :
- -117.0
- false_easting :
- 500000.0
- false_northing :
- 0.0
- scale_factor_at_central_meridian :
- 0.9996
- spatial_ref :
- PROJCS["UTM Zone 11, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
- GeoTransform :
- 699960.0 30.0 0.0 4100040.0 0.0 -30.0
array(0)
- _FillValue :
- -9999.0
- scale_factor :
- 0.0001
- add_offset :
- 0.0
- long_name :
- Red
da_red.hvplot.image(x='x', y='y', cmap='gray', aspect='equal')Unable to display output for mime type(s):